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EXAMINATION OF TWO METHODS OF DESCRIBING THE 
THERMODYNAMIC PROPERTIES OF OXYGEN 
NEAR THE CRITICAL POINT 

By Thomas H. Rees and John T. Suttles 
Langley Research Center 

SUMMARY 

A computer study was conducted to compare the numerical behavior of two 
approaches to describing the thermodynamic properties of oxygen near the critical point. 
One approach utilized properties from the tabulated data presented by L. A. Weber, and 
the other approach was based on computing properties using the curve-fit equations 
developed by R. B. Stewart. Data on the relative differences between values of specific 
heats at constant pressure (Cp^, density, and isotherm and isochor derivatives of the equa- 
tion of state are presented for selected supercritical pressures at temperatures in the 
range 100 to 300 K. The results of a more detailed study of the c p representations 
afforded by the two methods are also presented. 

The correlation of the two methods at the selected pressures was very good outside 
the near-critical temperature region (150 to 170 K). Large differences (up to 30 percent) 
occurred within this region, however. The detailed study of Cp representations illus- 
trated that these large discrepancies were due not to differences in the basic data of the 
methods but to the fact that the resolution of Weber's data in the near-critical region was 
not adequate to allow sufficient description of the thermodynamic functions by linear 
interpolation. In order to illustrate the effects of this poor description, identical calcu- 
lations were performed based on the two methods of obtaining thermodynamic properties. 
Regions of erratic behavior in the calculations based on Weber's tabulations were cor- 
related to errors introduced by linear interpolation. 

INTRODUCTION 

The removal of the mixing fans from the supercritical oxygen storage tanks in the 
Apollo spacecraft introduced the possibility that the fluid in the tank could be stratified, 
possibly giving rise to several serious effects. First, it could cause erroneous quantity 
probe measurements with the potential for precipitating a mission abort. Second, if the 
low-density fluid from the stratified layer entered the oxygen supply lines, the potential 
would exist for fuel-cell shutdown. Finally, if the stratified layer became too extensive, 



a perturbing force of even extremely low magnitude could cause depressurization with 
the potential for the existence of two -phase fluid in the tank. 

The concern about the effects of stratification in the redesigned tanks led to a sig- 
nificant volume of analytical study in that area (refs. 1,2, and 3). The thermodynamic 
data for oxygen that served as bases for the various studies were not identical; some 
investigators used the tabulated data presented by L. A. Weber (ref. 4), while some relied 
on the curve -fit equations developed by R. B. Stewart (ref. 5). 

Comparisons of results presented in references 1 to 3 with unpublished calcula- 
tions made by J. T. Suttles and G. L. Smith of the Langley Research Center have revealed 
significant differences. For example, Suttles and Smith using Weber's tabulations for the 
thermodynamic data predict the generation of very rapid oscillations of the mean pres- 
sure in the tank (about 25 cycles per hour after 2 hours) for a high-density case whereas 
Baldwin, Reinhardt, and Sheaffer (ref. 1, ch. 6) using Stewart's equations predict rates 
of about 3^- cycles per hour. Although the overall methods of analysis in the cited papers 
are diverse and thus undoubtedly contribute to the varying results, the contributions to 
these differences caused by the use of different descriptions of the thermodynamic data 
were also questioned. The purpose of this paper is, therefore, to compare the descrip- 
tions of the thermodynamic data for supercritical oxygen afforded by the two methods 
mentioned (refs. 4 and 5) and to compare the results of identical calculations based upon 
the two methods. 

The approach taken was first to compute by each method values of specific heats 
at constant pressure (Cp), density, and the isotherm and isochor derivatives of the equa- 
tion of state for supercritical pressures of 58.5 and 60 atm over the temperature range 
100 to 300 K. The percent differences of values computed by Stewart's equations from 
interpolated values of Weber's tabulations were computed and plotted against tempera- 
ture. The values of Cp were also plotted along selected isobars and isotherms. 

Finally, thermodynamic properties computed by Stewart's equations were applied to a 
simple mathematical model of the Apollo cryogenic oxygen storage tank for which the 
results of calculations based on Weber's tabulations were available. 

SYMBOLS 

Aj coefficients in vapor-pressure curve fit (eq. (4)), defined in appendix B, 

table IIIB 

Ci coefficients in zero-pressure specific-heat equation (eq. (B3)), defined in 

appendix B, table IIB 
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specific heat at constant pressure, j/kg-K 

specific heat at zero pressure, liter-atm/mole-K 

Cp as a function of T only, that is, along an isobar 

functions for evaluating derivatives of equation of state (eq. 1)), 
defined in appendix A, table IVA 

functions for evaluating integrals of thermodynamic equations, 
defined in appendix B, table IB 

enthalpy, liter -atm/mole 

reference enthalpy at triple point, liter -atm/mole 
mass flux, kg/hr 

coefficients in Stewart's curve fit to equation of state for oxygen (eq. (1)), 
defined in appendix A, table IA 

pressure, atm (1 atm = 101.325 kN/m2) 

vapor pressure, atm 

characteristic gas constant for oxygen, 0.0820535 liter-atm/mole-K 
temperature, K 

temperature at critical point, K 
temperature at triple point, K 
time, sec 

internal energy, liter -atm/mole 
heat flux, watts 
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functions for evaluating equation of state and derivatives, defined in 
appendix A 

functions for evaluating integrals of thermodynamic equations, defined 
in appendix B 

density, mole/liter 

critical density, mole/liter 

saturated liquid density, mole/liter 

saturated vapor density, mole/liter 
isochor derivative of equation of state 

isotherm derivative of equation of state 


Subscripts: 


COND 

conduction 

HTR 

heater 

input 

input 

leak 

leak 

net 

net 

RAD 

radiation 


Xi 

Yi 

P 


Pc 

p sat L 

p saty 



METHOD 

Two methods of describing the thermodynamic properties of oxygen were utilized 
in the present investigation. The first method was based on the least-squares curve fit 
to the equation of state for oxygen developed by Stewart (ref. 5). Stewart’s curve fit 
gives pressure in terms of temperature and density over the entire range of data avail - 
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able at the time. The estimated accuracies of the equation of state over the three main 
temperature -pressure ranges are given as follows: 

1. Gaseous range of values (65 to 300 K, to 340 atm) - ±0.1 percent in p 

2. Liquid range (65 K to T c , to 340 atm) - ±0.15 percent in p 

3. Near-critical range (125 to 165 K, 40 to 100 atm) - ±0.3 percent in p 

Stewart also presented methods for calculating entropy, enthalpy, internal energies, and 
specific heats from derived functions using his equation of state and its derivatives. 

The second method of describing the thermodynamic properties of oxygen utilized 
in the present investigation involved interpolation of- Weber's tabulations of the thermo- 
dynamic properties (ref. 4). Values of molar volume, isotherm and isochor derivatives 
of the equation of state, enthalpy, entropy, specific heats, and velocity of sound were 
tabulated as functions of pressure (to 330 atm) and temperature (T 0 to 300 K). These 
data were presented in increments of 5 atm and 2 K up to 160 K and in increments of 
5 atm and 5 K above 160 K. 

For the purpose of the present study, values of Cp, p, anc * ("ir) were 

taken directly from Weber’s tabulations (ref. 4), read into the digital computer, and ref- 
erenced by linear interpolation. 


The equations required f±r utilization of Stewart's method were adapted from ref- 
erence 5. The equations and procedures used in the present investigation are presented 
in detail. Also, a computer listing of subroutines used for applying Stewart's method is 
included as appendix C. 


Equation of State 


The curve fits for the thermodynamic properties of oxygen developed by Stewart 
(ref. 5) were programed into the computer. The least-squares curve fit for the equation 
of state, given in reference 5 as equation (7), is presented here as equation (1): 


p = pRT + 


n l T + n 2 



n g T z + n 


7 T + n 8 



+ ( n ll T + n l 2 )P + i n 13 + "t~J P + p 


n 14L5 . 3ri5 . "16 , n 17 


exp 


( n 25 p2 ) 


+ P 


5 




n 24 p 


n 28 +1 




exp 






(1) 
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The coefficients for equation (1), n^ through n2g, are included in table IA of appen- 
dix A. The constraints for the critical point and the fixed -point data imposed by Stewart 
upon the curve fit for the equation of state are listed in tables HA and IHA of appendix A, 
respectively. 

Since the calculations required density as a function of temperature and pressure, 
equation (1) was solved numerically for density by an iterative method. 

The isochor derivative and the isotherm derivative of equation (1) 

are given in appendix A as equations (A2) and (A3), respectively. 

Specific Heat Calculations 

Since the method of Stewart (ref. 5) does not include a procedure for calculating Cp 
throughout the entire range of supercritical conditions being considered, it was necessary 
to adopt the procedure described in this section. The approach is based upon the use of 
the enthalpy equation to produce a consistent calculation of Cp for the entire range of 
interest. 

The specific heat at constant pressure was approximated by a central differencing 
of the enthalpy versus temperature curve; 

H(T + AT, pi) - H(T - AT ,p 9 ) 

hr (2) 

where pj and p 2 are densities at T + AT and T - AT, respectively, for a given 
pressure. A temperature increment AT = 1.0 K was used in equation (2), except in the 
near-critical region (150 KIT = 170 K) where AT = 0.1 K was used to improve the 
approximation. A brief computer study showed that, for AT = 0.1 K, a halving of the 
temperature increment to AT = 0.05 K produced a change in c p of less than 0.2 percent 
in the near -critical region, which is within the accuracy of Stewart's curve -fit equations. 

Using the method of reference 5, it was necessary to calculate the enthalpies -for 
equation (2) differently for temperatures above and below T c . The procedure is 
described in the following sections. 

Temperatures above critical .- For temperatures above T c the enthalpy used in 
equation (2) was calculated by equation (18) of reference 5, given here as equation (3) 



6 



The solution to the integrals over p in equation (3) is simplified by using the equation 
of state (eq. (1)) and its isochor derivative (eq. (A2)). The indefinite solutions to the three 
integrals in equation (3) are given in appendix B. The reference enthalpy used was the 
enthalpy of the ideal gas at the triple point T 0 = 55.0 K which is = 15.7 liter- 

atm/mole (ref. 5). 

Temp eratures below critical . - The method used for calculation of the enthalpy of 
supercritical oxygen at temperatures below T c involved integration of internal energies 
from the saturation point along an isotherm. The use of internal energies in this region 
is part of the method of reference 5. The sequence of computations is described as 
follows: 


First, the vapor pressure was calculated from the vapor -pressure curve fit for 
oxygen given as equation (7) in reference 4, 

log e p y (T) = A 1 + A 2 T + A 3 T 2 + A 4 .T 3 + A 5 T 4 + AgT 5 + A ? T 6 + AgT 7 (4) 

The A^’s for equation (4) are listed in table IIIB of appendix B. 

Next, the saturated liquid and vapor densities, p ga j. and Psaty’ were calculated 

by the simultaneous solution of the equation of state (eq. (1)) and the vapor-pressure 
curve (eq. (4)). 


The saturated -vapor enthalpies were calculated by evaluating equation (3) at T, 
p ga t v > and Py The saturated -liquid enthalpies were calculated by subtracting the 

enthalpy change due to vaporization as given by the Clapeyron equation from the saturated- 
vapor enthalpy: 


(5) 


H(T,P satL ) = H ( T,p satv) -T(^)^-^ 

/dp v \ 

where 1 derivative of the vapor -pressure curve (eq. (4)). 

The internal energies and the enthalpies were related as follows: 

U(T,p) = H(T,p)--E 

The internal energies were calculated by equation (24) of reference 5, given here as 
follows: 


(6) 


U(T,p) = U( T ,p satL ) - T «jp + £ 


' p sat L P 2 


(7) 


7 



where u ( T >P sa t L ^ is the saturated -liquid internal energy, calculated from equations (5) 

and (6). The solutions to the integrals in equation (7) were facilitated by the use of the 
following relations: 



After the appropriate substitutions are made, equation (7) becomes 



The integrals used in equation (10) are the same as those used in equation (3). 

Finally, the enthalpy used in equation (2) was calculated by equation (6). 

Numerical difficulties . - Iterating for the saturation densities was complicated by 
the existence of an extraneous positive real root of equation (1) which lies between the 
roots representing the saturated vapor and liquid densities. Calculation of the extraneous 
root was avoided by careful choice of the bounds on the iteration. Evaluation of equa- 
tion (A3), the isotherm derivative of the equation of state, at each calculated root pro- 
vided a check for extraneous roots, since the isotherm derivative is always positive at 
the saturation densities but negative at the extraneous root. The problem of the extrane- 
ous root to equation (1) is addressed more comprehensively by Walter A. Reinhardt in 
reference 1, chapter 4. 

Discontinuities were calculated in the Cp -isobars at T c . Since the temperature 
range encompassing the discontinuities was only about 2 K at the pressures investigated, 
this problem was circumvented by calculating Cp's above and below the discontinuity 
and interpolating linearly along the isobar for intermediate temperatures. 

Stewart's equations (ref. 5) exhibited sensitivity to propagation of errors. The 
large number of operations involved tend to mushroom seemingly insignificant relative 
errors in p, for example, into large relative errors in p calculated from equation (1). 
This effect can be minimized by demanding higher -order accuracy on the iterations for p, 
and by programing for minimum propagation of errors. (Reduction of the number of com- 
puter operations is one effective method for minimizing propagated errors.) A discus- 
sion of computer propagation of errors may be found in a numerical -methods text such 
as reference 6. 
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The method described above for calculation of Cp's proved to be expensive in 
computer run time relative to a table look-up method; consequently, care was required 
in the coding and logical flow to insure maximum efficiency. 

PROCEDURE 

The procedure used in the present study involved first calculating values of Cp, p, 

, and (— ) for two supercritical pressures (58.5 and 60 atm) and a range of tem- 
T V^/p 

peratures, 100 to 300 K, in 2 K increments. The percent differences between values 
from Stewart's method and from Weber's tabulations were computed by the following 
relation: 

Percent difference = Stewart's. value. -.Weber' s value x 1Q0 

Weber's value 

These percent differences were plotted against temperature in figure 1. Values of Cp 
were then calculated over the temperature range 130 to 180 K in increments of 0.2 K 
along the isobars 55, 60, and 58.5 atm and are presented as plots in figure 2. In figure 3, 
a family of similar curves from reference 4 is presented for comparison. Values of c 

r 

were also calculated over a range of pressures, 55 to 70 atm, in increments of 1 atm for 
four isotherms, 140, 157, 160, and 180 K. These isotherms are presented as figure 4. 

In order to explore the effects on calculated results caused by using Stewart's 
method rather than Weber's, identical calculations were performed and compared. In 
these calculations use was made of a simple mathematical model of the Apollo cryogenic 
oxygen storage tank developed by J. T. Suttles and G. L. Smith. Suttles and Smith used 
a one -dimensional cylindrical tank model as shown in figure 5 which facilitated a solution 
but retained the essential features of the actual tank. Their calculations gave the time 
rate of change of pressure as a function of temperature for two mass flux rates, 

0.2268 kg/hr or 0.5 lb/hr and 0.6804 kg/hr or 1.5 lb/hr, both for heat addition (Q n et = 

125 watts) and for no heat addition. The present calculations were made for pressures 
of 58.5 and 60 atm. The 60-atm case shows comparison of the results based on the two 
thermodynamic representations at a pressure at which Weber's data were tabulated. The 
58.5-atm case duplicates the conditions imposed on the calculations of Suttles and Smith. 
The results of the calculations are presented as plots in figure 6. 

RESULTS AND DISCUSSION 

Thermodynamic Properties 

The percent differences in the thermodynamic properties calculated by the different 
methods (figs. 1(a) to 1(d)) range between zero and 3 percent outside the near -critical 

9 




region (150 K < T < 170 K); however, large differences were noted inside this region. 

Also, the error introduced by linear interpolation of Weber’s tabulations appears to have 
a significant effect upon the percent differences. The magnitudes of the percent differ- 
ences were, in general, larger at points representing interpolated Weber's data than at 
points representing tabulated data. The differences in Cp were by far the most signifi- 
cant; therefore, Cp was selected for a detailed examination. 

Figures 2(a) to 2(c) present comparisons of values of Cp produced by the differ- 
ent methods plotted as functions of temperature along three isobars (55 atm, 60 atm, and 
58.5 atm). It should again be noted that agreement outside the near-critical temperature 
region (150 K to 170 K) appears to be very good. Furthermore, the points corresponding 
to tabulated data in reference 4 (circled points, 55 and 60 atm) agree very well with values 
calculated by Stewart's curve fits over the entire temperature range; if the circled points 
were faired smoothly, the faired curves would closely approximate the curves calculated 
from Stewart's equations. This statement is supported by figure 3, which was taken 
directly from reference 4. Although the isobars from figure 2 do not appear in figure 3, 
the smoothness and characteristic shape of the 50- and 70-atm isobars of figure 3 give 
strong support to at least the shape of the Cp -curves obtained from Stewart's equations 
shown in figure 2. 

Since figures 2(a) and 2(b) represent data at pressures at which values of c p were 
tabulated in reference 4, the poor behavior of Weber's data in the proximity of the peaks 
could only be due to linear interpolation of Cp(T) between points too widely spaced in 
temperature to describe the peak adequately. Figure 2(c), on the other hand, represents 
the 58.5-atm isobar, which is between tabulated pressures. This figure, therefore, illus- 
trates the combined effect of double first-order interpolation of Cp, which appears to be 
only slightly worse than the effect of interpolating Cp(T) along the tabulated isobars. 

Figure 4, a comparison of Cp’s along selected isotherms, is presented to illus- 
trate better the effects of linear interpolation of pressure. The 140 and 180 K isotherms 
lie outside the near -critical region and reassert the statement that first-order interpola- 
tion of Weber's data is sufficient outside of this region. The 160 K isotherm (a tabulated 
temperature) illustrates the irregularity introduced by linear interpolation of Cp with 
pressure in regions where the slope of the curve changes rapidly. Although the relative 
accuracy of the two methods at points where tabulated data exist in reference 4 has not 
been established, certainly Stewart's equations yield smoother results in the near- 
critical region than linear interpolation of Weber's tabulated data. Finally, the 157 K 
isotherm again represents a comparison of the results of Stewart's method with double- 
interpolated Weber's data corresponding to the peak of figure 2(a). The greatest part of 
the deviation of this pair of curves near 55 atm is attributed to interpolation of Cp with 
temperature since 55 atm is a tabulated pressure. The error caused by first-order 
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interpolation of c p with pressure appears to be relatively insignificant compared to 
the gross error caused by first-order interpolation of the Cp(T) -curves near their peaks. 


Calculations for Tank With Uniform Properties 

The calculations for a uniform tank performed by Suttles and Smith were repeated, 
and the results are presented in figures 6(a) and 6(b). The calculations gave the time 
rate of change of pressure as a function of temperature for two mass flux rates, 

0.2268 kg/hr or 0.5 lb/hr and 0.6804 kg/hr or 1.5 lb/hr, both for heat addition ^Q n et = 

125 watts) and for no heat addition. The calculations were made using both methods for 
description of the thermodynamic properties of oxygen for two pressures (p = 60 atm 
and 58.5 atm). Results from the two methods agreed over the entire temperature range 
for no heat addition, and for temperatures greater than 170 K for the heat -addition case. 
The methods also agreed fairly well for the latter case up to approximately 140 K. The 
divergent behavior between 140 and 154 K correlates with differences in the thermody- 
namic properties previously noted in that range. 

The relative behavior of the results based upon the two methods in the region 
between 154 and 170 K warrants special attention, however. The results from Weber's 
data in figure 6(a) up to 160 K correspond to tabulated thermodynamic properties from 
reference 4. It is important to note that the pressure derivatives calculated from Weber's 
data at 162 and 164 K correspond to regions of poor description of the peak in the c p (T)- 
curve by linear interpolation (figs. 2(b) and 2(c)), while the derivatives at 58.5 atm, 156 
and 158 K, correspond to points of large percent differences in all the thermodynamic 
properties (fig. 1). 

In summary, comparison of figures 6(a) and 6(b) indicates that the variations in the 
pressure derivative between 154 and 170 K calculated from linearly interpolated data 
from reference 4 (these variations were described as "unusual" by Suttles and Smith) 
were caused by the error introduced by such interpolation. 

The calculations for the stratified -tank problem performed by Suttles and Smith 
were not made, but significant differences in the behavior of calculations around the criti- 
cal temperature could be expected. 

CONCLUDING REMARKS 

A study was made to compare the approximations of the thermodynamic behavior of 
supercritical oxygen as afforded by Weber's tabulations and by Stewart's curve -fit equa- 
tions as presented in this paper. Comparisons were presented for density, the isotherm 
and isochor derivatives of the equation of state, and specific heat at constant pressure 
(Cp) for a range of temperatures and two supercritical pressures. Special attention was 
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given to approximations of c p . Calculations for a simple one -dimensional mathematical 
model of the redesigned Apollo cryogenic oxygen storage tank with uniform properties ' 
were made based upon both approximations of the thermodynamic properties of oxygen. 
The results were plotted and presented for comparison. 

In general, the thermodynamic properties obtained by the two methods agreed very 
well except in the near -critical temperature region (150 to 170 K). Stewart’s equations 
were more effective in approximating the general shapes of the curves of the thermody- 
namic properties in regions of rapid variation than linear interpolation of Weber's tabu- 
lations. The resolution of Weber's data was insufficient for first-order interpolation in 
regions of rapid variation in the thermodynamic properties. Additional tabulated data 
near the critical point would greatly enhance the reliability of Weber's tabulations in this 
region. 

Similarly, the results of the uniform -tank calculations agreed for temperatures 
above 170 K. Regions of erratic numerical behavior in the results of the calculations 
based on Weber's tabulations were correlated to errors introduced by linear interpolation. 

From a utilization standpoint, both methods of describing the thermodynamic prop- 
erties of oxygen have advantages. For limited ranges of temperature and pressures such 
as those in f he present study, the table look-up method was considered more convenient 
and efficient for computer programing. On the other hand, if a very wide range of con- 
ditions must be studied, the large number of data points required for this method could 
make computer storage requirements the limiting factor; the curve -fit method becomes 
more advantageous as the ranges of temperature and pressure increase. 

If Stewart's equations are to be employed, certain obvious techniques can be incor- 
porated to enhance the efficiency of the method; for instance, the problem should be for- 
mulated so that temperature and density are the independent variables. Owing to the 
large number of operations involved, calculation of specific heats should be avoided if 
possible, especially in the near-critical region where the uncertainty of the representa- 
tion is the greatest. Stewart's equations should be programed to avoid redundant opera- 
tions and to minimize propagation of errors. 

More recent data are now available (refs. 7 and 8) which improve the description 
of the thermodynamic behavior of oxygen in the critical region; nevertheless, linear inter- 
polation in this region is not recommended. Regardless of the methods of describing 
thermodynamic properties chosen for future investigations, a brief numerical study of the 
description afforded will probably disclose any gross misrepresentations of the data 
similar to linear interpolation of Weber's tabulations in the critical region. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., July 25, 1972. 
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APPENDIX A 


DERIVATIVES OF THE EQUATION OF STATE (REF. 5) 


Equation of State 


The equation of state (eq. (1)) may be written 


where the 

24 

p = pRT + ^ njXj 
i=l 

n^'s are the coefficients listed in table LA, 

and the Xj functions are as 

follows: 

x x = p 2 T 

X 9 = p 3 /t 

x }7 = p3 fl /T 4 


X 2 = P 2 

x 10 = p 3 /t 2 

Xl8 = P 5 «l/T 2 


x 3 = P 2/t 2 

x u = p 4 T 

x 19 - P V t3 


X 4 = p 2 / T 4 

Xj 2 = P 4 

x 20 = P 5 J l/T 4 


x 5 = p 2 /t 6 

X i3 = P 5 

x 2 i = p 7 «i/t 2 


X 6 = p 3 T 2 

X i4 = P 5 / T 

x 22 = P 1f l/T 3 


x 7 = p 3 T 

x 15 = P 3l l/ T2 

x 23 = P 7f l/T 4 


x 8 = p 3 

x 16 = p \/ t 3 

X 2 4=P (n28+ \f 3 

The fj's 

are given in table IVA. 



Isochor Derivative 


The isochor derivative for equation (Al) is given as 

\9T; p 

24 

d)„ - pR + I 


(A2) 


where the n^’s are the coefficients listed in table LA, and the Xj functions are as 
follows: 


13 



APPENDIX A - Continued 


Xj = p 2 

X 9 = -p 3 /t 2 

X 17 = - 4 P V T5 

x 2 = 0 

X 10 = -2p 3 /t 3 

x 18 = - 2 P 5 fl/T 3 

x 3 = - 2 p 2 /t 3 

X 11 = p4 

x 19 = -3p 5 fl/T 4 

X 4 = -4p 2 /T 5 

Xi2 = 0 

X 20 = -> 5 fl/T 5 

X 5 = - 6 P 2 /t 7 

Xi 3 = 0 

X 21 = Vfl/T 3 

X 6 = 2Tp 3 

x 14 = - p 5 /t 2 

x 22 = - 3 P 7 fl/T 4 

X 7 = p 3 

x 15 - Vfl/ T 3 

X 23 = -4P VT 5 

x 8 = o 

Xl 6 = -3p 3 fi/T 4 

x _ n ( n 28 +1 ) f f 
x 24 “ P f 2 f 4 


The fj’s are given in table IVA. 


Isotherm Derivative 


The isotherm derivative ( — ] for equation (Al) is given by 

Wt 


V 

d pj T 


24 

= RT + ^ npq 
i=l 


(A3) 


where the n> are the coefficients listed in table IA, and the Xj functions are as 
follows: 


Xj = 2pT 

X 9 = 3p 2 /t 

x n = ' 6 / t4 

X 2 = 2 p 

X 1Q = 3p2/T 2 

X 18 = f?/T 2 

X 3 = 2p/T 2 

X n = 4p 3 T 

X19 - f 7 / T3 

X 4 = 2p/t 4 

X u = 4p 3 

x 20 = '7/t 4 

X 5 = 2 P /t 6 

X 13 = 5p 4 

X21 = %/t 2 

X 6 = 3p 2 T 2 

x 14 = 5 P 4 / T 

x 22 = %/ t3 

X 7 = 3p 2 T 

X 15 = f 6 /T 2 

x 23 = %/ t4 

X 8 = 3p 2 

x i6 = f eA 3 

x 24 = f 9 f 2 f 3 


>28+1) ( n 28 +1 ) 
> I10I3 + P 


f 2 f ll 


The f> are given in table IVA. 
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APPENDIX a - Continued 


TABLE IA.- COEFFICIENTS FOR EQUATION (1) 


T 

in K, p in atm, p in mole/liter 



R = 0.0820535 

n 1Q = -3.59419602 X 10 

n 19 = 

-2.67817667 x 10 2 

n x = 3.38759078 x 10' 3 

n n = 1.02209557 X 10" 6 

n 20 = 

1.05670904 x 10 5 

n 2 = -1.31606223 

n 12 = 1.90454505 X 10' 4 

n 21 = 

5.63771075 X 10" 3 

n 3 = -7.38828523 X 10 3 

n 13 = 1.21708394 x 10" 5 

n 22 “ 

-1.12012813 

n 4 = 1.92049067 X 10 7 

n 14 = 2.44255945 x 10' 3 

n 23 = 

1.46829491 x 10 2 

n 5 = -2.90260005 x 10 10 

n 15 = 1.73655508 x 10 2 

n 24 = 

9.98868924 x 10~ 4 

n 6 = -5.70101162 X 10~ 8 

n 16 = 3.01752841 x 10 5 

n 25 = 

-0.00560 

n ? = 7.96822375 x 10' 5 

n 1? = -3.49528517 x 10 7 

n 26 = 

-0.157 

n 8 = 6.07022502 X 10‘ 3 

n 18 = 8.86724004 x 10' 1 

n 27 = 

-0.350 

n 9 = -2.71019658 


n 28 = 

0.90 


TABLE HA.- CONSTRAINTS IMPOSED ON EQUATION (1) 
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APPENDIX A - Concluded 


TABLE HIA.- FIXED -POINT DATA 


Critical pressure 

50.14 atm 

Critical temperature 

154.77 K 

Critical density 

13.333 mole/liter 

Normal boiling temperature (International practical 
temperature scale (IPTS), fixed point) 

90.18 K 

Density saturated vapor at normal boiling point (eq. (1)) 

0.1396 mole/liter 

Density saturated liquid at normal boiling point (eq. (1)) 

35.65 mole/liter 

Triple point pressure 

0.00150 atm 

Triple point temperature 

54.353 K 


TABLE IVA.- FUNCTIONS FOR EVALUATING DERIVATIVES 
OF EQUATION OF STATE 


f 1 = exp(n 25 p 2 ) 

f 10 = n 28 p ^ 28 ^ 

f 2 = p"28 - Pc "28 

fll = 2f 2 n 26 f 3 f 10 

f 3 = exp[n 26 f 2 2 + n 27 (T - T c ) 2 ] 

f 12 = 2f 5 pn 25 + 2f l n 25 

f 4 = 2n 27 (T - T c )f 3 

f 1 3 = 6f lP + 6f 5 p 2 + f 12P 3 

f 5 = 2f lP n 25 

f 14 = 20fj P 3 + 10f 5P 4 + f 12 /3 5 

f 6 = 3f lP 2 + f 5P 3 

f 15 = 42f lP 5 + 14f sP 6 + f 12P 7 

f 7 = 5f lP 4 + fg/) 5 

f 16 = n 28 ( n 28 + J ) p 

% = 7f lP 6 + fjp 7 

f - n (n - l)o^ n28 “ 2 ^ 
t 17 _n 28l n 28 1 ) p 

f 9 = ( n 28 + X ) p 28 

f 18 = ( 2f 10 n 26 f 3 + 2f 2 n 26 f ll) f 10 + 2f 2 n 26 f 3 f 17 
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APPENDIX B 


INTEGRATION OF THE THERMODYNAMIC EQUATIONS (REF. 5) 

The solution of the integral ^ ^ dp (from eq. (3)) at constant tempera- 
ture is given by equation (Bl). This integration uses the isochor derivative given 

\ 8T /p 

in appendix A as equation (A2). 

r- _ 24 

ffe-AfS Op=y n.Y, (Bl) 


__ 

( S. . J-fiJL'll dp = y n Y 

J[P p2WpJ T 4 


where the n^'s are the coefficients listed in table IA, and the Yj functions are as 
follows: 


Yi= -P 

Yg = p 2 /(2T 2 ) 

y 17 = 4g 4 /T 5 

o 

n 

CM 

>1 

Y 10 - P 2 A 3 

y 18 = W t3 

Y 3 = 2p/T 3 

Y u = -p3/3 

Y 19 - 3 8 5 / t4 

&■ 

ii 

h" 

U1 

Yl2 = 0 

y 20 = 4 «sA 5 

Y 5 = 6p/t7 

y 13 = ° 

y 21 " 2g 6 / T3 

y 6 = p 2 t 

y 14 - P 4 /(«T 2 ) 

y 22 - 3B 6i /T 4 

Y 7 = -p 2 / 2 

y 15 = ^/r 3 

y 23 = 4 S 6 / t5 

>< 

00 

1! 

o 

y 16 = 3g 4 /T« 

Y 24 = " n 27( T 


where the g^s are given in table IB. 


The solution of the integral 


i h-% 


jdp at constant temperature is given by 


\K f 

equation (B2). This integration uses p = f(T,p) from equation of state (eq. (Al)) as given 
in appendix A. 

24 

ife - dp - L n ‘ Y ‘ (B2) 

\P / T i=l 
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APPENDIX B - Continued 


where the n^’s are the coefficients listed in table LA, and the Y^ functions are as 
follows: 

Y 1 = Tp 

Y 9 = p 2 /(2T) 

Y 17 = g 4 A 4 

Y 2 = P 

Y 10 = P 2 /(2T 2 ) 

Y 18 = SsA 2 

y 3 = p/t 2 

Y u =p 3 T/3 

Y i9 - e 5 A 3 

y 4 = p/t 4 

Y 12 = p 3 /3 

Y 20 - SsA 4 

Y 5 = p 2 T6 

Yi3 = p 4 /4 

Y 21 = g 6 A 2 

Y 6 = p 2 T 2 /2 

Y h = p 4 /(4T) 

Y 2 2 = g 6 /T 3 

Y 7 = p 2 t/2 

Y 15 = g4/ T2 

Y 23 = g 6 / T4 

Q. 

11 

oo 

Y 16 = S 4 / t3 

Y 24 = g 3/( 2n 28 n 26) 


where the g^'s are given in table IB. 

The solution of the integral \ c° dT is given by 

tr 

8 

J c£ dT = R £ (B3) 

i=l 

where R is given in table IA, and the C^’s are the coefficients listed in table IIB, and 
the Yj functions are as follows: 

Yj = -l/(2T 2 ) Y 5 = T 2 /2 

Y 2 = -1/T Y 6 = T 3 /3 

Y 3 = log e T Y 7 = T 4 /4 

Y 4 = T Yg = uT^e 11 - l) 

where u = C g /r. 
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APPENDIX B - Concluded 


TABLE IB. - FUNCTIONS FOR EVALUATING INTEGRALS 
OF THERMODYNAMIC EQUATIONS 



TABLE IIB.- COEFFICIENTS FOR ZERO-PRESSURE SPECIFIC-HEAT EQUATION 


C 1 = -1.86442361 X 10 2 

Cg = -1.11035799 X 10' 8 

C 2 = 2.07840241 x 10 

C 7 = 2.08612876 x 10' 11 

Cg = -3.42642911 X 10' 1 

Cg = 1.01894691 

C 4 = 3.50297163 

C 9 = 2.23918105 x 10 3 

C 5 = 2.05866482 X 10’ 7 



TABLE IIIB.- COEFFICIENTS FOR VAPOR -PRESSURE EQUATION (EQ. 4) 

T in K, Py in atm, p in mole/liter 
Aj = -62.5967185 
A 2 = 2.47450429 
A 3 = -4.68973315 x 10' 2 
A 4 = 5.48202337 x 10" 4 
Ag = -4.09349868 x 10 -6 
A 6 = 1.91471914 x 10‘ 8 
A 7 = -5.13113688 x 10“ U 
A g = 6.02656934 X 10" 14 
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APPENDIX C 


LISTING OF FORTRAN SUBROUTINES FOR STEWART’S METHOD 


SUBROUTINE STE WART < TEMP . RHO , PRESS , CP * DPDT , DPDRHO . ENERGY ) 

C INPUT PRESS (ATM) RHO (GM-MOLE/L ) 

DIMENSION XN<28) 

C0MMON/N/XN,TEMPCR,RHOCRIT*R 

DATA (XN( I) .1 = 1-. 28 ) /3. 387590 78E-3 « -1 . 3 1 606223 , -7 . 38828523E+3 , 
1 1 .92049067E+7,-2.90260005E+10,-5.70101 1 62E-8 , 7 . 96822375E-5 , 
26. 07022502E-3.-2.71019658.-3.59419602E+I » I.02209557E-6, 
31.90454505E-4, l.2l708394E-5',2.44255945E-3, 1.736555 08E + 2* 
43.0175284] E + 5, -3. 495285 17E* 7,8. 86724004E-1 , -2 . 678 1 7667E+2 , 

51 .05670904E + 5. 5.63771 075E-3.-1 . 12012813, 1 .46829491 E + 2* 

69. 98868924E-4,-. 00560,-. 157, -.350, .9/ 

RHODUM=RHO 

PCR I T = 50 .14 

01 =C 1 C2C4 ( TEMP , RHO) 

DPDT=C2< TEMP, RHO) *1. 01325E*5 
DPDRH0=C4 (TEMP, RHO) * 1 . 0 1 325E+5/3 1 . 9988 
CP=CPHU (TEMP, RHO, 01 ) * 1 . 0 1 325E + 5/31 . 9988 

PRESS=01*1.01 325E+5 
ENERGY=0 . 

RHO=RHOOUM 

RETURN 

END 
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APPENDIX C — Continued 


FUNCTION CPHU(TEMP.RHO.P) 

ALL PRESSURES IN ATM ALL DENSITIES IN GM-MOLE/L 
CPHU (L-ATM/GM-MOLE/K) H ( L-A TM/GM-MOLE ) U.DITTO 
EXTERNAL FOE X 
COMMON/THERM/PI .T 
COMMON C 1 «C2 
DIMENSION A ( 8 ) .H (2) 

DATA <A ( I ) , 1 = 1 , 8 ) /6.0265693AE-14.-5. 131 136R8E-1 1 * 1 . 9 147 1 9 1 AE-8 . 

1- 4. 09349868E-6. 5. 48202337E-4, -4. 689733 15E-2. 2. 47450429, 

2- 62.5967185/ 

ENTHALP (T.RHO.RHOl .P) = 1590. 929*9. 86896E -3+ P/RHO-.0820535»T 
A-05 (55. »RHO) 

1 -D2 ( T , RHOl ) -D1D2D5 < T » RHOl ) *T ♦D5 ( T , RHO ) *D2 (T.RHO) ♦D1D2D5(T,RH0)*T 
ENERGY (T.RHO.RHOl .USAT) =USAT-D2 ( T, RHO 1 ) -D1D2D5( T.RHO 1 ) *T+D2< T.RHO) 
1+D1D2D5 (T.RHO) *T 

R=0. 0820535 % RHODUM=RHO $ ICOUNT=0 $EREL=5.E-7 $ TINCR=1.0 
T I = 1 53 . S TIP1=155.2 

IFtTEMP.GT.TI. AND. TEMP. LT.TIPDGO TO 2222 
GO TO 1111 
2222 FACTOR=31. 9988/1. 01325E5 
CP I =C 1 *F ACTOR 
CPIP1 =C2*FACTOR 

CP=CPI* (TEMP-TI)/(TIP1-TI) *( CPI PI -CPI) 

CPHU=CP 

RETURN 

11 1 1 IF (TEMP.GE. 150. . AND.TEMP.LE. 170. ) TINCR=0. 1 
DELT=2 . *T I NCR 

6 IF (ICOUNT.EQ.2)GO TO 205 

T=TEMP-T INCR % ICOUNT=ICOUNT*l $ TINCR=-TINCR 
ASSIGN 1 TO INDEX 
NDEX= 1 
P 1 — P 

RGM=.8*RH0$RGP=1.2*RH0*DELTRG=. 1*RGM 
IF (T.GT. 155. .AND.T.LT. 162. ) RGM=.5*RH0 
EREL=5 . E- 7 
EA 8 S=ERFL' tt RGM 

9 CALL ITR2 (RHO.RGM.RGP.DELTRG.FOFX.EREL. EABS ♦ 50 . I CODE ) 

IF ( I CODE . NE . 0 ) GO TO 1000 
GO TO INDEX. (1.55.60) 

1 IF (T.LT.154.77)G0 TO 10 
5 H ( ICOUNT) =ENTHALP (T.RHO. 0. ,P) 

GO TO 6 
10 PVLN=A(1) 

DO 15 1=2.8 
15 PVLN-PVLN*T+A ( I ) 

PV=EXP (PVLN) 

IF (P.LE.PV) GO TO 5 
RS AVE = RHO $ P1=PV 

59 pgm=A.«(EXP(0.0232*PV)-1.) $ RGP=RH0*.6$ DELTRG=. 1*RGP-. 1*RGM 
IF (T.GT. 153. ) EREL=5. E- 1 0 

EA 8 S=EPEL»RGM 

NDEX=60 

ASSIGN 60 TO INDEX 
GO TO 9 

60 RSATV=RHO 

RGM=30.-1.5*RSATV $ RGP=1 . 1*RSAVE S DELTRG= ( RGP-RGM ) * . 2 

EABS=EREL*RGM 

NDEX=55 



APPENDIX C - Continued 


ASSIGN 55 TO INDEX 
GO TO 9 
55 RSAT L = RHO 
90 CONTINUE 
SUM=A(i)*7. 

DO 100 1=2*7 
100 SUM=SUM*T+ (S-I ) *A ( I ) 

DPVD T = SUM*PV 

DELTH=T*OPVDT* ( 1 . /RSAT V- 1 . /RSA TL ) 

150 HPSATV=ENTHALP(T*PSATV.O. . PV) 

hrsatl=hrsatv-oelth 

usat=hrsatl-pv/rsatl 

RHO=RS AVE 

U=ENERGY (T fRHO.RSATL.USAT) 

H ( I COUNT ) = U+-P/RHO 
GO TO 6 

205 CPHU= (H (2) -H ( l ) ) /DELT 

TI = TEMP-TINCR % T2 = T E M P ♦ T I NC R 

900 FORMAT ( 10X ,*T1 = »*F8.3.* T2 = * ,F8.3»» H(T1) = **E11.4* 

1' * H(T2)=**E11.4»* RV*RL *»2E11.4) 

RHO=RHOOUM 

RETURN 

1000 WRITE (6*2) ICODE*NDEX 

^ FORMAT ( 1H0 *21HERR0R IN CPHU IC0DE=I2*5X**INDEX = »*I2///) 
return 
END 
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APPENDIX C - Continued 


FUNCTION C1C2C4(TEMP*RH0) 

DIMENSION X ( 24 ) * XN ( 28 ) .F (18) 

COMMON/N/XN* TEMPCR . RHOCR I T . R 
R=0. 0820535 
RHOCR I T= 1 3 . 333 
TEMPCR= 154.77 

F(l)=EXP(XN(25)*RHO»'>2) S F (2) =RHO**XN (28 ) -RHOCRI T**XN (28) 

F (3) =EXP(XN(26) *F (2) **2+XN (27) »(TEMP-TEMPCR)**2) 

X ( 1 ) =RH0**2»TEMP $ X(2)=RhO»*2 $ X (3) =RH0**2/TEMP**2 

X (4) =RH0**2/TEMP**4 $ X (5) = RH0**2/TEMP*«-6 $ X (6) =RHO**3*TEMP*»2 

X ( 7) =RHO**3»TEMP $ X(8)=RHO**3 S X <9 ) =RH0**3/TEMP 

X ( 10) =RHO**3/TEMP»*2 $ X ( 11 ) =RH0**4*TEMP $ X(12)=RH0»*4 

X ( 1 3 ) =RHO**5 S X(14) =RhO**5/TEMP $ X ( 1 5 ) =RH0*"3*F ( 1 ) /TEMP**2 

X(16), = RH0*»3»F(1)/TEMP**3 $ X ( 1 7 ) =RH0**3*F ( 1 ) /TEMP**4 

X ( 18) =RH0**5*F ( 1 ) /TEMP**2 $ X ( 19) =RH0**5*F ( 1 ) /TEMP**3 

X (20) =RHO**5*F (1) /TEMP»<w* $ X ( 2 1 ) =RH0**7*F ( 1 ) /TF.MP«*2 

X (22) =RHO**7*F ( 1 ) /TEMP**3 $ X (23) =RHO*»7*F ( 1 ) /TEMP**4 

X(24)=RH0**(XN(28) ♦ 1 . ) *F (2) *F (3) 

FACTOR=RHO*R*TEMP 
GO TO 10 
ENTRY C2 

F (4) =2.*XN(27) MTEMP-TEMPCR) *F (3) 5 F (5) =2. *F ( 1 ) »RHO»XN (25) 
X(l)=RHO*»2 $ X ( 2 ) =0 . % X ( 3) =-2. »RH0**2/TEMP**3 
X (4) =-4.*PH0**2/TEMP»*5 $ X (5) =-6. «RH0**2/TEMP»*7 

X (6) =2. <> TEMP*RHO**3 S X(7)=RH0**3 % X(8)=0. $ X ( 9 ) =-RHO**3/TEMP**2 
X ( 10) =-2.»RHO**3/TEMP**3 $ X ( 1 1 ) =RH0**4 $ X < 1 2) =0 . $ X < 1 3 ) = 0 . 

X ( 14) =-RH0**5/TEMP**2 $ X ( 15) =-2 . *RH0**3*F < 1 ) /TEMP**3 
X ( 16) =-3.*RHO**3*F ( 1) /TEMP**4 $ X ( 1 7) =-4. *RH0**3 tt F ( 1 ) /TEMP**5 
X ( 18) =-2.*RH0 <m>5*F ( 1) /TEMP#*3 $ X ( 1 9 ) =-3 . *RH0**5*F ( 1 ) /TEMP**4 
X (20 ) =-4.«RH0**5*F ( 1 ) /TEMP-»*5 $ X ( 21 ) =-2 . *RHO** 7*F ( 1 ) /TEMP**3 
X (22) =-3. *RHO**7*F ( 1 ) /TEMP**4 $ X (23) =-4. »RHO**7*F ( 1 ) /TEMP**5 
X (24) =RH0«* (XN (28) ♦ 1 . ) »F (2) *F (4) 

FACTOR=RHO*R 
GO TO 10 
ENTRY C4 

F (6) =3.*F ( 1) »RH0**2*F (5) *RHO**3 S F (7) =5.*F ( 1 ) *RH0**4*F (5) *RH0*»5 
F (8) =7. *F ( 1) *RH0*»6*F (5) *RHO**7 $ F (9) = (XN (28) ♦ 1 . > *RHO**XN (28) 
F(10)=XN(28)*RHO**(XN(28)-1.) $ F < 1 1 ) =2 . *F ( 2 ) «XN < 26 ) *F < 3) *F < 1 0 ) 

X ( 1 ) =2.*RH0*TEMP * X ( 2 ) =2. *RHO $ X ( 3) =2 . *RH0/TEMP*«2 
X(4)=2.»RH0/TEMP**4 S X (5) =2.*RH0/TEMP**6 $ X ( 6 ) =3 . *RHO**2*TEMP**2 
X(7)=3.*RHO**2*TE^P S X (8) =3. *RHO**2 $ X (9) =3.#RHO**2/TEMP 
X ( 10) =3.*RH0**2/TEMP»*2 $ X ( 1 1 ) =4 . *RH0*»3»TEMP i X ( 1 2 ) =4 . »RHO**3 
X ( 13) =5.*RHO»*4 $ X(14)=5.*RH0**4/T£MP $ X ( 15) =F (6) /TEMP**2 
X ( 16) =F (6) /TEMP**3 S X ( 17) =F (6) /TEMP«*4 $ X ( 1 8 > =F < 7 ) /TE MP**2 
X < 19) =F (7) /TEMP**3 $ X (20) =F ( 7) /TEMP**4' S X (21 ) =F (8) /TEMP**2 
X (22) =F (8) /TEMP**3 $ X (23) =F (8) /TEMP*»4 

X (24) =F (9) *F (2) »F (3) ♦RHO** ( XN ( 28 ) ♦ 1 . ) *F ( 1 0 ) "F ( 3 ) +KH0** ( XN ( 28 ) +1 . ) * 
IF (2) *F (11) 

FACTOR=R*TEMP 
10 SUN=0. 

DO 1 1=1*24 
1 SUM=SUM+XN(I)*X(I) 

C1C2C4=SUM+FACT0R 

RETURN 

END 
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APPENDIX C — Concluded 


FUNCTION D 1 D2D5 ( T , RmO) 

DIMENSION XN (28) »G (6) »Y (24) , C (9) 

COMMON /N/XN ,TCR,RHOCR,R 

DATA (C ( I ) • 1=1 ,9) /-1 .86442361E + 2, 2. 07840241E+1 *-3.4264291 IE-1 , 

1 3.50297 1 6.3 <2. 05366482E-7 , - 1 . 1 1 035799E-8 , 2. 086 1 2876E- 1 1 * 
21.01894691 ,2. 2391 81 G5E+3/ 

T E M P C R = 7 C R $ Ri-IOCR I T =RHOCRi>TEMP=T 

G(1)=EXP(XN(25)*RH0»*2) * G < 2 ) =RH0**XN ( 28 ) -RHOCRI T**XN ( 28 ) 
G(4)=G(1)/(2.*XN(25) ) 

G (5) =G ( 1 ) *(XN(25) *RH0»*2-1.)/ ( 2 . *XN ( 25) **2 ) 

G (6) =G( 1 ) * (XN(25) **2*RH0**4-2.* (XN(25) *RH0**2-1 . ) ) / ( 2 . *XN ( 25 ) ) 

G (3) =E XP (XN(26) #G (2) »»2+XN ( 27) # ( TEMP-TEMPCR) **2> 

RHOSQ=RHO*RHO * T2 = T*T S T3 = T2*T $ T4=T3»T $ T5=T4»T 
Y(l)=-RHO $ Y < 2 ) =0 . $ Y (3) =2.“RH0/T3 $ Y (4) =4.*RH0/T5 

Y (5) =6.*RHO/T**7 S Y(6)=RH0SU*T i Y ( 7 ) =-RH0SQ/2 . $ Y(8)=0. 

Y (9)=RhOSG/(2. <t T2) $ Y { 1 0 ) =RH0SQ/T3 $ Y < 1 1 ) =-RH0*.*3/3 . 

Y ( 12) =Y ( 13) =0. S y(14)=RH0**4/(4.*T2> S Y ( 15) =2.*G (4) /T3 
Y ( 16) = 3.*G (4) /T4 $ Y (1 7 ) =4 . *6 (4 ) / T5 S Y ( 1 8) =2 . *G ( 5 ) /T3 
Y(19)=3.*G(5)/T4 S Y<20)=4.*G(5)/T5 $ Y ( 21 ) =2. *G ( 6) /T3 

Y (22) = 3 . *G ( 6 ) /T4 <5 Y ( 23 ) =4 . »G ( 6 ) /T5 

Y(24) = -XM27) 4(T-TCR)4G(3)/(XN(28)4XN(26) > 

GO TO 10 
ENTRY 02 

PH03=RH0SQ*RH0 $ RH04=RH03*RHO 

Y(1)=T4RH0 % Y (2) =RH0 S Y(3)=RH0/T2 S Y(4)=RH0/T4 $ Y ( 5) =RHO/T **6 
Y (6) =RH0SQ»T2/2. % Y < 7 ) =RHOSu4T/2 . $ Y ( 8 ) =RHOSQ/2 . 

Y (9) =.5»RH0SQ/T $ Y(10)=.5*RHOSQ/T2 $ Y ( 1 1 > =RH03*T/3 . 

Y ( 12) =RH03/3. <5 YU3)=RH04/4. $ Y ( 14) = . 25*RH04/T 
Y(15)=G(4)/T2 % Y(16)=G(4)/T3 $ Y ( 17) =G (4) /T4 
Y(18)=G(5)/T2 $ Y ( 19) =G ( 5 ) /T 3 $ Y (20) =G (5) /T4 

Y (21 ) =G (6) /T2 $ Y (22) =G(6) /T3 S> Y ( 23) =G ( 6) /T4 
Y(24)=G(3)/(2. <I ’XN(28)4XN(26) ) 

10 SUM=0. 

DO 11 1=1,24 

11 SUM=SUM*Y(I)"XN(I) 

01 0205= SUM 
return 

ENTRY 05 

U=C(9)/T S T2=T tt T $ T3=T2*T S> T4=T3*T 

Y ( 1 ) =- . 5/T2 S Y(2)=-l./T i Y ( 3 ) = ALOG ( T ) $ Y(4)=T $ Y(5)=.5*T2 
Y (6 ) =T 3/3 . $ Y (7) =T4/4 . S Y (8) =U*T/ (EXP ( U> -1 . ) 

S U M = 0 . 

DO 20 1=1,8 
20 SUM = SUM + C ( I ) *Y ( I ) 

010205= SUM *R 

RETURN 

END 
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(a) Specific heat c p . 

Figure 1.- Percent difference between values of thermodynamic properties 
calculated by Stewart's method and values from interpolation (Weber's 
tabulations). Percent differences given by 

Stewart's value - Weber's value x 
Weber's value 
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(b) Density p. 
Figure 1.- Continued. 
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(c) Isochor derivative 



Figure 1.- Continued. 
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(d) Isotherm derivative 



Figure 1.- Concluded. 
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(a) p = 55 atm. 

Figure 2.- Comparison of specific heat at constant pressure along selected isobars. 
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(a) p = 60 atm, a pressure tabulated by Weber. 

Figure 6.- Comparison of uniform -tank calculations using Stewart’s and Weber's thermodynamic properties. 
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